load moviedata;

M = 3952;
N = 6040;
D = 30;

lambda = 0.01;
epsilon = 0.001;
maxepoch=80;
 
mean_rating = mean(train_vec(:,3)); 
rating = train_vec(:,3);
rating = rating - mean_rating;
rating = rating';
nr = length(rating);

u_ind = train_vec(:,1);
i_ind = train_vec(:,2);

U = 0.1*randn(D,N);
V = 0.1*randn(D,M);

ts = tic;
for epoch=1:maxepoch
    pred_train = sum(U(:,u_ind).*V(:,i_ind),1);
    train_rmse = sqrt(mean((pred_train-rating).^2));
   
    pred_probe = sum(U(:,probe_vec(:,1)).*V(:,probe_vec(:,2)),1);
    probe_rmse = sqrt(mean((pred_probe+mean_rating-probe_vec(:,3)').^2));
   
    fprintf('epoch %3d %2.4f %2.4f\n',epoch,train_rmse,probe_rmse);  
    
    I0 = repmat(rating - sum(U(:,u_ind).*V(:,i_ind),1),D,1);
    IV = I0.*V(:,i_ind) ;
    IU = I0.*U(:,u_ind) ;
   
    dU = zeros(D,N);
    dV = zeros(D,M);
    
    for ii = 1:length(rating)
       dU(:,u_ind(ii)) =  dU(:,u_ind(ii)) + IV(:,ii);
       dV(:,i_ind(ii)) =  dV(:,i_ind(ii)) + IU(:,ii);
    end
    
    dU = dU - lambda*U;
    dV = dV - lambda*V;
    
    U = U + epsilon*dU;
    V = V + epsilon*dV;
end
toc(ts)
